import numpy as np
from scipy import interpolate, integrate
import matplotlib as mpl
import matplotlib.pyplot as plt


def plot(title=None):
    mpl.rcParams['font.sans-serif'] = ['SimHei']  # 指定默认字体
    mpl.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题
    if title:
        plt.title(title)
    plt.show()


def np_interp():
    """
    numpy 插值
    :return:
    """
    x = np.linspace(0, 8 * np.pi, 20)
    y = np.sin(x / 2)

    x_vals = np.linspace(0, 8 * np.pi, 100)
    # 一维线性插值
    y_interp = np.interp(x=x_vals, xp=x, fp=y)
    # 同上,kind为插值类型
    # 阶梯插值（0阶B样条曲线）：zero, nearest
    # 线性插值(1阶)：slinear, linear
    # 多阶B样条曲线：quadratic, cubic
    interpolate.interp1d(x, y, kind='linear')

    plt.plot(x, y, 'o')
    plt.plot(x_vals, y_interp, '-x')

    plot(title='numpy 插值')


def sp_interp():
    """
    numpy scipy 插值
    :return:
    """
    t = np.arange(0, 2.5, 0.1)
    x = np.sin(2 * np.pi * t)
    y = np.cos(2 * np.pi * t)

    # 计算欲插值曲线的样条系数
    # 参数s确定平滑点数，通常是m-sqrt(m),其中m是曲线点数；s=0表示不需要平滑
    # 输出第一个数：曲线点；第二个数：计算出来的系数；
    tcktuples, uarray = interpolate.splprep([x, y], s=0)
    unew = np.arange(0, 1.01, 0.01)
    # 进行样条插值计算（B样条和它的导数进行插值）
    splinevalues = interpolate.splev(unew, tcktuples)

    plt.figure(figsize=(10, 10))
    #  原始点
    # plt.plot(x, y, 'x', x, y, 'b')
    #  三次样条插值
    # plt.plot(splinevalues[0], splinevalues[1])
    #  真实数据
    # plt.plot(np.sin(2 * np.pi * unew), np.cos(2 * np.pi * unew))

    #  三图一块画
    plt.plot(x, y, 'x', x, y, 'b',
             splinevalues[0], splinevalues[1],
             np.sin(2 * np.pi * unew), np.cos(2 * np.pi * unew))

    plt.legend(['线条', '三次样条', 'True'])
    plot(title='numpy scipy 插值')


def sp_trape():
    """
    积分计算方法
    对f(x)在[a, b]上求积分，本质上：求x=a, x=b, y=0, y=f(x)围成曲边梯形的面积，在双低存在的情况下，计算高即可
    难度：梯形变为曲线

    1、梯形规则(trapezoidal rule)：
        (f(a) + f(b)) * (b - a) / 2

    2、辛普森规则(Simpson rule)：
        (f(a) + 4 * f((b - a) / 2) + f(b)) * (b - a) / 6

    3、复合梯形规则(composite trapezoidal rule)：
        (f(a) + sum(f(x(i)) * (b - a) / 2) + f(b)) * (b - a) / (2 * n)  # i为0到n-1

    4、复合辛普森规则(composite simpson rule)：
        (f(a) + 4 * sum(f(x((i + 1) / 2))) + 2 * sum(f(x(i))) + f(b)) * (b - a) / (6 * n)  # i同上

    5、龙贝格求积公式(Romberg)：迭代计算

    scipy 数值积分
    本例：计算 9-x^2 的从-3 到 3 的积分
    :return:
    """
    def source_func(a):
        return 9 - np.power(a, 2)

    def trape_func(a):
        # 积分
        return 9 * a - np.power(a, 3) / 3

    begin = -3
    end = 3
    N = 10

    x = np.linspace(begin, end, N)
    y = source_func(a=x)
    yromb = lambda x: source_func(a=x)

    # 三种方法计算积分值
    # np.trapz() : 使用复合梯形法则沿给定轴进行积分值计算
    """
    np.trapz([1,2,3]) = 4.0  
    计算点[1, 1], [2, 2], [3, 3]三个点与x轴围成的梯形面积
    (1 + 3) * (3 - 1) / 2 = 4
    
    np.trapz([1,2,3], x=[4, 6, 8]) = 8.0    
    计算点[4, 1], [6, 2], [8, 3]三个点与x轴围成的梯形面积
    (1 + 3) * (8 - 4) / 2 = 8
    
    np.trapz([1,2,3], dx=2) = 8.0    
    计算点[1, 1], [3, 2], [5, 3]三个点与x轴围成的梯形面积
    (1 + 3) * (5 - 1) / 2 = 8
    
    """
    t = np.trapz(y, x)
    # 同np.trapz
    t_sp = integrate.trapz(y, x)

    # composite Simpson’s rule
    s = integrate.simps(y, x)

    r = integrate.romberg(yromb, begin, end)

    # 计算真实实际的积分值
    aiv = trape_func(a=end) - trape_func(a=begin)

    info = 'trapz: {}\nsimps: {}\nromberg: {}\naiv: {}'.format(t, s, r, aiv)
    print(info)


sp_trape()
